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Abstract Many interesting and fundamentally practical optimization prob- 
lems, ranging from optics, to signal processing, to radar and acoustics, involve 
constraints on the Fourier transform of a function. It is well-known that the fast 
Fourier transform (fft) is a recursive algorithm that can dramatically improve 
the efficiency for computing the discrete Fourier transform. However, because 
it is recursive, it is difficult to embed into a linear optimization problem. In 
this paper, we explain the main idea behind the fast Fourier transform and 
show how to adapt it in such a manner as to make it encodable as constraints 
in an optimization problem. We demonstrate a real-world problem from the 
field of high-contrast imaging. On this problem, dramatic improvements are 
translated to an ability to solve problems with a much finer grid of discrctized 
points. As we shall show, in general, the "fast Fourier" version of the optimiza- 
tion constraints produces a larger but sparser constraint matrix and therefore 
one can think of the fast Fourier transform as a method of sparsifying the 
constraints in an optimization problem, which is usually a good thing. 
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1 Fourier Transforms in Engineering 

Many problems in engineering involve maximizing (or minimizing) a linear 
functional of an unknown real-valued design function / subject to constraints 
on its Fourier transform / at certain points in transform space ([I]). Examples 
include antenna array synthesis (see, e.g., p~2 | fT3 | rT6 ] ) . FIR filter design (see, 

e.g., HCHM]), and corona g ra P h desi § n ( see > e -s-> IsifTWnjUPiirn^nirTiigifn] ) . 

If the design function / can be constrained to vanish outside a compact interval 
C = (—a, a) of the real line centered at the origin, then we can write the Fourier 
transform as 

7(0 = P e 2 "^f(x)dx 



and an optimization problem might look like 
maximize J" c(x)f(x)dx 

subject to — e < < e, £ e D , , 

-e<9f/(0<e, ^ D 
< f(x) < 1, xeC, 

where D is a given subset of the real line, e is a given constant, and 3?(z) and 
9(2) denote the real and imaginary parts of the complex number z. In Section 
[7J we will discuss a specific real-world problem that fits a two-dimensional 
version of this optimization paradigm and for which dramatic computational 
improvements can be made. 

Problem is linear but it is infinite dimensional. The first step to making 
a tractable problem is to discretize both sets C and D so that the continuous 
Fourier transform can be approximated by a discrete Riemann sum: 

n 

% = e MkAx ^f k Ax, -n<j< n. (2) 

k — — n 

Here, n denotes the level of discretization, 

2a 

Ax 



2n + 1 



A£ denotes the discretization spacing in transform space, /& = f(kAx), and 

Computing the discrete approximation ([2| by simply summing the terms 
in its definition requires on the order of ./V 2 operations, where N = 2n + 1 is 
the number of discrete points in both the function space and the transform 
space (later we will generalize to allow a different number of points in the 
discretization of C and D). 

Choosing A£ too large creates redundancy in the discrete approximation 
due to periodicity of the complex exponential function and hence one generally 
chooses At; such that 

AxAt < 1. 
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In many real-world applications, A£ is chosen so that this inequality is an 
equality: A£ = \/{NAx). In this case, the Riemann sum approximation is 
called the discrete Fourier transform. 

2 A Fast Fourier Transform 

Over the past half century there has been an explosion of research into algo- 
rithms for efficiently computing Fourier transforms. Any algorithm that can 
do the job in a constant times NlogN multiplications/additions is called a 
fast Fourier transform (see, e.g., [51l2"l ll5l l()]). There are several algorithms that 
can be called fast Fourier transforms. Here, we present one that applies natu- 
rally to Fourier transforms expressed as in ([2]). In this section, we assume that 
A^ = l/(NAx). 

A sum from — n to n has an odd number of terms: N = In + 1. Suppose, 
for this section, that N is a power of three: 

N = 3 m . 

Fast Fourier transform algorithms assume that it is possible to factor N into 
a product 

N = N Ni. 
For the algorithm of this section, we put 

iV = 3, and N 1 = 3 m ~ 1 . 

The first key idea in fast Fourier transform algorithms is to write the single 
sum Q as a double sum and simultaneously to represent the discrete set of 
transform values as a two-dimensional array of values rather than as a one- 
dimensional vector. Specifically, we decompose k as 

k = 7V A:i + k 

so that 

— n < k < n —no < ko < uq and — ni < k\ < rii, 

where 

n = {N - l)/2 = (3 - l)/2 = 1 

and 

n x = (N, - l)/2 = (3™- 1 - l)/2. 
Similarly, we decompose j as 

3 = Niji + j 

so that 

—n < j < n — n± < jo < n\ and — 1 < j\ < 1. 
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With these notations, we rewrite the Fourier transform |2| as a double sum: 

1 Tlx 

fjo,n = E E e^ N » k ^ A < N ^ + ^f koM Ax, (3) 



kg — — 1 ki— — rt\ 



where / fco , fel = /jv fci+fc and /jWi = fmji+jo ■ Distributing the multiplications 
over the sums, we can rewrite the exponential as 



e 2ni(N k 1 +k )Ax(N 1 j 1 +jo)A£ 



_ e 2-KiN k 1 Ax(N 1 j 1 +j )A£ e 2-Kik a Ax(N 1 j 1 +j Q )A(, 

_ e 2mN k 1 AxN 1 j 1 A£ e 2-KiN kiAxj A£ e 2TTik Ax(N l j l +] a )A£, 

_ ^TviNakxAxjaAi e 2mk Ax (N 1 j 1 +j )A£ 



where the last equality follows from our assumption that NqN\AxA^ = NAxAt; 
1. Substituting into ([3]), we get 

1 / ni \ 

j, j = E e 2 ' KikoAx( ' Nljl+: > ' >Ai E e 2mNoklAx3oA ^ ■ fk k Ax 

fco=— 1 \fcl=— til / 

We can compute this nested sum in two steps: 



9*M= E e 2 " N " klAxjoAi f koM Ax, "i|i° ^T' (4) 

fcl=— Tlx 
1 



f. . - \" p 2mk AxjAi 



ko=-l 



-«i < jo < «i, 

-i<h < 1. 



By design, computing /j ^ for —ni < j < n% and — 1 < ji < 1 is equivalent 
to computing /j for — n < j < n. 

2.1 Complexity 

If we compute fj 0> j x in two steps according to the equations given above, then 
the number of multiply/adds is 

N 2 N + NN = N(N X + N ). 

On the other hand, the one-step algorithm given by {2} requires N 2 multi- 
ply/adds. Hence, the two-step algorithm beats the one-step algorithm by a 
factor of 

N 2 N 

n(n 1 + n o) = n7Tn- ~ n ^ = n ° = 3 - 
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2.2 Recursive Application 

One can do better by iterating the above two-step algorithm. From the formula 
for gj ,k given in Q, we see that g is a discrete Fourier transform of a subset 
of the elements of the vector {fy : k — —n, . . . ,n} obtained by sampling 
/ at a cadence of one every Nq elements. And, the coefficient NqAxA^ in 
the exponential equals N /N = 1/Ni, which again matches the number of 
terms in the sum. Hence, we can apply the two-step algorithm again to this 
Fourier transform. The second key component of the fast Fourier transform is 
the observation that this process can be repeated until the Fourier transform 
only involves a sum consisting of a single term. 

Let //v denote the number of multiply/adds needed using the recursive 
algorithm to solve a problem of size N — 3 m . Keeping in mind that N = 3, 
we get 

I N = I 3m =3/ 3m -! +3-3™ 

= 3(3/ 3 ™-a + 3-3 m - 1 ) + 3 m+1 
= 3 2 / 3m -2 + 2 • 3 m+1 

= 3*J 3 m-» + k ■ 3 m+1 

= 3 m 7 3 o + m • 3 m+1 

= 3 m (l + 3m) 

= iV(l + 31og 3 N). 

Hence, the recursive variant of the algorithm takes on the order of iVlog 3 N 
operations. 

3 A General Factor-Based Algorithm 

The advantage of fast Fourier transforms, such as the one presented in the 
previous section, is that they have order TV log N complexity. But, they have 
disadvantages too. One disadvantage is the need to apply the basic two-step 
computation recursively. Recursion is fine for computing a Fourier transform, 
but our aim is to encode a Fourier transform within an optimization model. 
In such a context, it is far better to use a non-recursive algorithm. 

A simple modification to the two-step process described in the previous sec- 
tion produces a variant of the two-step algorithm that makes a more substan- 
tial improvement in the initial two-step computation than what we obtained 
before. The idea is to factor N into a pair of factors with each factor close to 
the square-root of N rather than into 3 and N/3. Indeed, in this section, we 
assume, as before, that N can be factored into 

N = N N X 
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but we do not assume that Nq — 3. In fact, we prefer to have Nq N\. As 
before, we assume that N — 2n + 1 is odd and therefore that both No and Ni 
are odd: 

N = 2n + l and Ni = 2m + 1. 

At the same time, we will now assume that the number of points in the 
discretization of the Fourier transform does not necessarily match the num- 
ber of points in the discretization of the function itself. In many real-world 
examples, the "resolution" of the one discretization does not need to match 
the other and artificially enforcing such a match invariably results in a slower 
algorithm. So, suppose that the discrete Fourier transform has the form 



k——n 



e 2nikAxjA tf k Ax, -m<j< m, (5) 



and let M = 2m+l denote the number of elements in the discretized transform. 
Again, M is odd and therefore we factor it into a product M — MqM\ of two 
odd factors: 

M = 2m + 1 and M\ = 2m 1 + 1. 
If we now decompose our sequencing indices k and j into 

k = N kx + k and j = M ji + jo, 

we get 



fjo,j 



e 2iriN k 1 AxMojiA£, ^iriNoktAxjoA^ e 2iTikoAx{M j 1 +j Q )A(, 



- E E 

ko =—no fei=-ni 

■fkoM Ax - 

As before, we need to assume that the first exponential factor evaluates to 
one. To make that happen, we assume that NqMoAxA^ is an integer. In real- 
world problems, there is generally substantial freedom in the choice of each of 
these four factors and therefore guaranteeing that the product is an integer is 
generally not a restriction. With that first exponential factor out of the way, 
we can again write down a two-step algorithm 

ni 

_ \ ^ 2-niNo^AxjaAZ r a ~m < Jo < m , 

9jo,k - 2^ e JkoM^' - n „<k< n „. 



k\ — — ni 
n 



no < k < no, 

2irik Q Ax(M j 1 +jo)A( n . , ~ m < jo < n%Q 

-mi < h < m i- 



? - \ " p 2irikaAx(Maj 1 +j„)A£ 

feo=— n 
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3.1 Complexity 

The number of multiply/adds required for this two-step algorithm is 
NM + MN = MN ( — + —). 

If M f» iV and Mi pa N\ ps viV, the complexity simplifies to 

2iV%/7V. 

Compared to the one-step algorithm, which takes N 2 multiply /adds, this two- 
step algorithm gives an improvement of a factor of yN /2. This first-iteration 
improvement is much better than the factor of 3 improvement from the first 
iteration of the recursive algorithm of the previous section. Also, if M is much 
smaller than AT, we get further improvement over the full N x N case. 

Of course, if Mo, Mi, No, and Ni can be further factored, then this two- 
step algorithm can be extended in the same manner as was employed for the 
algorithm of the previous section successively factoring M and N until it is 
reduced to prime factors. But, our eventual aim in this paper is to embed these 
algorithms into an optimization algorithm and so we will focus our attention 
in this paper just on two-step algorithms and not their recursive application. 

4 Fourier Transforms in 2D 

Many real-world optimization problems, and in particular the one to be dis- 
cussed in Section [7J involve Fourier transforms in more than one dimension. 
It turns out that the core idea in the algorithms presented above, replacing a 
one-step computation with a two-step equivalent, presents itself in this higher- 
dimensional context as well |18| . 

Consider a two-dimensional Fourier transform 

f(Z,v) = J J e 2 ^ +m) f{x,y)dydx 

and its discrete approximation 

n n 

fh,h= E E e 2 ^ x ^+y^)f kiM AyAx, -m<ji,j 2 <m, 

k\— — n k%=— n 

where 





= kAx, 


— n < k < n, 


Vk 


= kAy, 


—n < k < n, 




= m, 


—m < j < m, 


m 


= j &n, 


—m < j < m, 


fki ,k 2 


= f(xk!,yk 2 ), 


—n < ki, fc 2 < n 


f 71 ,32 




-m < ji,j 2 < m 
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Performing the calculation in the obvious way requires M 2 N 2 complex 
additions and a similar number of multiplies. However, we can factor the ex- 
ponential into the product of two exponentials and break the process into two 
steps: 



n 

9nM = Yl e 2 ™"^ f kl ,k 2 Ax, -m < j 1 <m,-n<k 2 < n, 

ki——n 
n 

fh,h= E ^ myk ^ 2 9nM^y, -m<j u j 2 <m, 

k 2 = — n 



It is clear that, in this context, the two-step approach is simply to break up 
the two-dimensional integral into a nested pair of one-dimensional integrals. 
Formulated this way, the calculation requires only MN 2 + M 2 N complex 
additions (and a similar number of multiplications). 

The real-world example we shall discuss shortly involves a two-dimensional 
Fourier transform. Given that the idea behind speeding up a one-dimensional 
Fourier transform is to reformulate it as a two-dimensional transform and 
then applying the two-step speed up trick of the two-dimensional transform, 
we shall for the rest of the paper restrict our attention to problems that are 
two dimensional. 



5 Exploiting Symmetry 

Before discussing real-world examples and associated computational results, 
it is helpful to make one more simplifying assumption. If we assume that / is 
invariant under reflection about both the x and y axes, i.e., f(—x, y) = f(x, y) 
and f(x,—y) — f(x,y) for all x and y, then the transform has this same 
symmetry and is in fact real- valued. In this case, it is simpler to use an even 
number of grid-points (N = 2n and M = 2m) rather than an odd number and 
write the straightforward algorithm for the two-dimensional discrete Fourier 
transform as 



n n 

fh,h = 4 E E cos ^ x k 1 (j 1 )cos(2TTy k . 2 r] j2 )f klM AyAx 1 1 < j u j 2 < m, 
fe 1= i fc 2 =i 

(6) 
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where 



^ K 


= (k- 1/2) Ax, 


1 < k < n, 


Vk 


= (k-l/2)Ay, 


1 < k < n, 




= (j - 1/2)A£, 


1 < j < m, 


Vj 


= (j - 1/2)A V , 


1 < j < m, 


fki ,k2 


= f(xk 1 ,Vk 2 ), 


1 < k\ , fc 2 < n 


fjl J2 




1 < juh < m - 



The two-step algorithm then takes the following form: 

n 

9jiM = 2 C0S ( 2nx ki^i)fk 1 M Ax ^ 1 < ii <m,l<k 2 <n, 

k 1 = l 
n 

fjiA = 2 cos(2TTy k2 r} j2 )g jl!k2 Ay, 1 < j u j 2 < m, 
fe 2 =i 

5.1 Complexity 

The complexity of the straightforward one-step algorithm is m 2 n 2 and the 
complexity of the two-step algorithm is mn 2 + m 2 n. Since m = M/2 and 
n = N/2, we see that by exploiting symmetry the straightforward algorithm 
gets speeded up by a factor of 16 and the two-step algorithm gets speeded 
up by a factor of 8. But, the improvement is better than that as complex 
arithmetic has also been replaced by real arithmetic. One complex add is the 
same as two real adds and one complex multiply is equivalent to four real 
multiplies and two real adds. Hence, complex arithmetic is about four times 
more computationally expensive than real arithmetic. 

6 Matrix Notation 

As Fourier transforms are linear operators it is instructive to express our al- 
gorithms in matrix/vector notation. In this section, we shall do this for the 
two-dimensional Fourier transform. To this end, let F denote the nxn matrix 
with elements /fe ll / C2 , let G denote the m x n matrix with elements gj ly k 2 , let 
F denote the m x m matrix with elements fj 1 j 2 , and let K denote the mx n 
Fourier kernel matrix whose elements are 

K ji,k 2 = cos(2wx kl ^j 1 )Ax. 



For notational simplicity, assume that the discretization in y is the same as it 
is in x, i.e., Ax = Ay, and that the discretization in rj is the same as it is in £, 
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i.e., Ar) = A£. Then, the two-dimensional Fourier transform F can be written 
simply as 

F = KFK T 

and the computation of the transform in two steps is just the statement that 
the two matrix multiplications can, and should, be done separately: 

G = KF 

F — GK T . 

When linear expressions are passed to a linear programming code, the 
variables are passed as a vector and the constraints are expressed in terms of 
a matrix of coefficients times this vector. The matrix F above represents the 
variables in the optimization problem. If we let fk, k = 1, . . . , n denote the n 
columns of this matrix, i.e., F = [fx f 2 ■ ■ ■ /„], then we can list the elements 
in column- by-column order to make a column vector (of length n 2 ): 



vec(F) 



fi 
h 

In 



Similarly, we can list the elements of G and F in column vectors too: 



vec(G) 



.9i 

92 



and 



It is straightforward to check that 



vec(G) 



K 



K 



and that 



vec(F) 



K 2,ll K 2,2l 



vec(F) 



fi 

h 



K 



vec(F) 



vec(G), 



where I denotes an m x m identity matrix. 

The matrices in these two formulae are sparse: the first is block diago- 
nal and the second is built from identity matrices. Passing the constraints to 
a solver as these two sets of constraints introduces new variables and more 
constraints, but the constraints are very sparse. Alternatively, if we were to 



Fast Fourier Optimization 



11 



express vec(F) directly in terms of vec(F), these two sparse matrices would 
be multiplied together and a dense coefficient matrix would be passed to the 
solver. It is often the case that optimization problems expressed in terms of 
sparse matrices solve much faster than equivalent formulations involving dense 
matrices even when the latter involves fewer variables and/or constraints (see, 
e.g., EI]). 

7 A Real- World Example: High-Contrast Imaging 

Given the large number of planets discovered over the past decade by so- 
called "indirect" detection methods, there is great interest in building a special 
purpose telescope capable of imaging a very faint planet very close to its much 
brighter host star. This is a problem in high-contrast imaging. It is made 
difficult by the fact that light is a wave and therefore point sources, like the 
star and the much fainter planet, produce not just single points of light in the 
image but rather diffraction patterns — most of the light lands where ray-optics 
suggests it will but some of the light lands nearby but not exactly at this point. 
In a conventional telescope, the "wings" of the diffraction pattern produced by 
the star are many orders of magnitude brighter than any planet would be at 
the place where the planet might be. Hence, the starlight outshines the planet 
and makes the planet impossible to detect. But, it is possible to customize the 
diffraction pattern by designing an appropriate filter, or a mask, to put on the 
front of the telescope. While it is impossible to concentrate all of the starlight 
at the central point — to do so would violate the uncertainty principle — it is 
possible to control it in such a way that there is a very dark patch very close 
to the central spot. 

Suppose that we place a filter over the opening of a telescope with the 
property that the transmissivity of the filter varies from place to place over 
the surface of the filter. Let f(x, y) denote the transmissivity at location (x, y) 
on the surface of the filter ((0, 0) denotes the center of the filter). It turns out 
that the electromagnetic field in the image plane of such a telescope associated 
with a single point on-axis source (the star) is proportional to the Fourier 
transform of the filter function /. Choosing units in such a way that the 
telescope's opening has a diameter of one, the Fourier transform can be written 
as ^ ^ ^ ^ 

m,V)= [ [ e 2 ^ + y^f(x,y)dydx. (7) 

J -1/2 J -1/2 

The intensity of the light in the image is proportional to the magnitude squared 
of the electromagnetic field. 

Assuming that the underlying telescope has a circular opening of diameter 
one, we impose the following constraint on the function /: 

f(x,y)=0 for x 2 + y 2 > (1/2) 2 . 

As often happens in real-world problems, there are multiple competing 
goals. We wish to maximize the amount of light that passes through the filter 
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and at the same time minimize the amount of light that lands within a dark 
zone T> of the image plane. If too much light lands in the dark zone, the 
telescope will fail to detect the planets it is designed to find. Hence, this latter 
objective is usually formulated as a constraint. This leads to the following 
optimization problem: 



maximize 



subject to 



// f(x,y)dydx 
^ 2 



fit,v) <£, (£,»?) €Z>, (8) 
f(x,y) =0, x 2 +y 2 > (1/2) 2 , 
0< f(x,y) < 1, for all x,y. 

Here, e is a small positive constant representing the maximum level of bright- 
ness of the starlight in the dark zone. Without imposing further symmetry 
constraints on the function /, the Fourier transform / is complex valued. 
Hence this optimization problem has a linear objective function and both 
linear constraints and convex quadratic inequality constraints. Hence, a dis- 
cretized version can be solved (to a global optimum) using, say, interior-point 
methods. 

Assuming that the filter can be symmetric with respect to reflection about 
both axes (in real- world examples, this is often — but not always — possible; see 
[3j for several examples), the Fourier transform can be written as 

,1/2 ,1/2 

=4/ / cos(2TTx£)coa(2Tryri)f(x,y)dydx. 
Jo Jo 

In this case, the Fourier transform is real and so the convex quadratic inequal- 
ity constraint in ^ can be replaced with a pair of inequalities, 

-V£</te,T»)<vs, 

making the problem an infinite dimensional linear programming problem. 

Figure [T] shows an AMPL model formulation of this problem expressed in 
the straightforward one-step manner. Figure [2] on the other hand, shows an 
AMPL model for the same problem but with the Fourier transform expressed 
as a pair of transforms — the so-called two-step process. 

As Figures[3j[4j and[5]show, the optimal solution for the two models are, of 
course, essentially the same except for the improved resolution in the two-step 
version provided by a larger value for n (n — 1000 vs. n = 150). Using LOQO 
[22] as the interior-point method to solve the problems, both versions solve in 
a few hours on a modern computer. It is possible to solve even larger instances, 
say n = 2000, if one is willing to wait a day or so for a solution. Ultimately, 
higher resolution is actually important because manufacturing these masks in- 
volves replacing the pixellated mask with a spline-fitted smooth approximation 
and it is important to get this approximation correct. 

Table [T] summarizes problem statistics for the two versions of the model as 
well as a few other size choices. Table [2] summarizes solution statistics for these 
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param pi := 4*atan(l) ; 
param rhoO := 4; 
param rhol := 20; 

param n := 150; # discretization parameter 

param dx := l/(2*n); 
param dy := dx; 

set Xs := setof {j in 0.5..n-0.5 by 1} j/(2*n); 
set Ys := setof {j in 0.5..n-0.5 by 1} j/(2*n); 
set Pupil := setof {x in Xs, y in Ys: x~2+y~2 < 0.25} (x,y); 

var f {x in Xs, y in Ys : x"2 + y~2 < 0.25> >= 0, <= 1, := 0.5; 

param m := 35; # discretization parameter 

set Xis := setof {j in 0. .m> j*rhol/m; 

set Etas := setof {j in 0..m} j*rhol/m; 

set DarkHole := setof {xi in Xis, eta in Etas: 

xi~2+eta~2>=rho0~2 kk 

xi~2+eta~2<=rhol~2 kk 

eta <= xi } (xi,eta); 

var fhat {xi in Xis, eta in Etas} = 

4*sum {(x,y) in Pupil} f [x,y] *cos (2*pi*x*xi) *cos (2*pi*y*eta) *dx*dy ; 

var area = sum {(x,y) in Pupil} f [x,y] *dx*dy; 

maximize throughput: area; 

subject to sidelobe_pos {Cxi, eta) in DarkHole}: fhat[xi,eta] <= 10" (-5) *f hat [0, 0] ; 
subject to sidelobe_neg {Cxi, eta) in DarkHole}: -10" (-5) *f hat [0 ,0] <= fhat [xi , eta] ; 

solve ; 



Fig. 1 AM PL model for discretized version of problem |8jl assuming that the mask is sym- 
metric about the x and y axes. The dark zone T> is a pair of sectors of an annulus with inner 
radius 4 and outer radius 20. The optimal solution is shown in Figure [3] 

Table 1 Comparison between a few sizes of the one-step model shown in Figure^ and a 
few sizes of the two-step model shown in Figure [2] The column labeled nonzeros reports 
the number of nonzeros in the constraint matrix of the linear programming problem and the 
column arith. ops. The One-Step-250x35 problem is too large to solve by LOQO, which is 
compiled for a 32-bit architecture operating system. 



Model 


n 


m 


constraints 


variables 


nonzeros 


arith. ops. 


One step 


150 


35 


976 


17,672 


17,247,872 


17,196,541,336 


One step 


250 


35 


* 




* 




Two step 


150 


35 


7,672 


24,368 


839,240 


3,972,909,664 


Two step 


500 


35 


20,272 


215,660 


7,738,352 


11,854,305,444 


Two step 


1000 


35 


38,272 


822,715 


29,610,332 


23,532,807,719 



same problems. These problems were run as a single thread on a GNU/Linux 
(Red Hat Enterprise Linux Server release 5.7) x86_64 server with dual Xeon 
X5460s cpus (3.16 GHz with 4 cores each), 32 GB of RAM and a 6.1 MB 
cache. 

Real telescopes have opennings that are generally not just open unob- 
structed disks but, rather, typically have central obstructions supported by 
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param pi := 4*atan(l) ; 
param rhoO := 4; 
param rhol := 20; 

param n := 1000; # discretization parameter 

param dx := l/(2*n); 
param dy := dx; 

set Xs := setof {j in 0.5..n-0.5 by 1} j/(2*n); 
set Ys := setof {j in 0.5..n-0.5 by 1} j/(2*n); 
set Pupil := setof {x in Xs, y in Ys: x~2+y~2 < 0.25} (x,y); 

var f {x in Xs, y in Ys : x"2 + y~2 < 0.25> >= 0, <= 1, := 0.5; 

param m := 35; # discretization parameter 

set Xis := setof {j in 0. .m} j*rhol/m; 

set Etas := setof {j in 0..m} j*rhol/m; 

set DarkHole := setof {xi in Xis, eta in Etas: 

xi"2+eta~2>=rho0"2 kk 

xi"2+eta~2<=rhol"2 kk 

eta <= xi } (xi,eta); 

var g {xi in Xis, y in Ys}; 

var fhat -Cxi in Xis, eta in Etas}; 

var area = sum {(x,y) in Pupil} f [x,y] *dx*dy; 

maximize throughput: area; 

subject to g_def {xi in Xis, y in Ys} : 

g[xi,y] = 2*sum {x in Xs : (x,y) in Pupil} 
f [x,y] *cos (2*pi*x*xi) *dx; 

subject to fhat_def {xi in Xis, eta in Etas}: 
fhat [xi, eta] = 2*sum {y in Ys} 
g [xi ,y] *cos (2*pi*y*eta) *dy ; 

subject to sidelobe_pos {(xi,eta) in DarkHole}: fhat[xi,eta] <= 10" (-5) *f hat [0 ,0] ; 
subject to sidelobe_neg {(xi,eta) in DarkHole}: -10" (-5) *f hat [0 ,0] <= fhat [xi , eta] ; 

solve ; 

Fig. 2 AM PL model reformulated to exploit the two-step algorithm. The optimal solution 
is shown in Figure [1] 

Table 2 Hardware-specific performance comparison data. The results shown here were 
obtained using the default value for all of LOQO's tunable parameters. It is possible to 
reduce the iteration counts to about 100 or less on all the problems by increasing the value 
of the epsdiag parameter to about le-9. 



Model 


n 


m 


iterations 


primal objective 


dual objective 


cpu time (sec) 


One step 


150 


35 


54 


0.05374227247 


0.05374228041 


1380 


One step 


250 


35 


* 


* 


* 


* 


Two step 


150 


35 


185 


0.05374233071 


0.05374236091 


1064 


Two step 


500 


35 


187 


0.05395622255 


0.05395623990 


4922 


Two step 


1000 


35 


444 


0.05394366337 


0.05394369256 


26060 



Fast Fourier Optimization 



15 








-20 -15 -10 -5 5 10 15 20 



Fig. 4 The optimal filter from the two-step model shown in Figurej^Jand a logarithmic plot 
of the star's image. 



spiders. It is easy to extend the ideas presented here to handle such situations; 
see [3]. 

As explained in earlier sections, the two-step algorithm applied to a one- 
dimensional Fourier transform effectively makes a two-dimensional represen- 
tation of the problem and applies the same two-step algorithm that we have 
used for two-dimensional Fourier transforms. It is natural, therefore to consider 
whether we can get more efficiency gains by applying the two-step algorithm 
to each of the iterated one-dimensional Fourier transforms that make up the 
two-step algorithm for the two-dimensional Fourier transform. We leave such 
investigations for future work. 
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Fig. 5 Close up of the two masks to compare resolution. 




Fig. 6 Logarithmic stretches are useful but can be misleading. Left: the image of the star 
shown in a linear stretch. Right: the same image shown in a logarithmic stretch. 
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